#include "stdafx.h"
#include "cv.h"
#include "cxcore.h"
#include "highgui.h"
#include "math.h"

void RGB2YUV(IplImage *src, IplImage *dst)
{
IplImage* img_y = cvCreateImage (cvGetSize (src),IPL_DEPTH_64F, 1);
IplImage* img_u = cvCreateImage (cvGetSize (src),IPL_DEPTH_64F, 1);
IplImage* img_v = cvCreateImage (cvGetSize (src),IPL_DEPTH_64F, 1);
CvScalar s;
double y,u,v;
for(int i=0;i<src->height;i++)
    for(int j=0;j<src->width;j++)
	{
		s=cvGet2D(src,i,j);
		y=0.257*s.val[2]+0.504*s.val[1]+0.098*s.val[0]+16;
		if (y>235)
			y=235;
		u=-0.148*s.val[2]-0.291*s.val[1]+0.439*s.val[0]+128;
		if (u>240)
			u=240;
		v=0.439*s.val[2]-0.368*s.val[1]-0.071*s.val[0]+128;
		if (v>240)
			v=240;
		s.val[0]=y;
		s.val[1]=u;
		s.val[2]=v;
		cvSet2D(dst,i,j,s);
	}
cvSplit(dst,img_y,img_u,img_v,0);
cvSaveImage("./TestImg/out_y.bmp",img_y);
cvSaveImage("./TestImg/out_u.bmp",img_u);
cvSaveImage("./TestImg/out_v.bmp",img_v);
}

void YUV2RGB(IplImage *src1,IplImage *src2,IplImage *src3, IplImage *dst)
{
IplImage* img_rgb = cvCreateImage (cvGetSize (src1),IPL_DEPTH_64F, 3);
CvScalar s,s1,s2,s3;
double r,g,b;
for(int i=0;i<src1->height;i++)
    for(int j=0;j<src1->width;j++)
	{
		s1=cvGet2D(src1,i,j);
		s2=cvGet2D(src2,i,j);
		s3=cvGet2D(src3,i,j);
		b=1.164*(s1.val[0]-16) + 2.018*(s2.val[0]- 128);
		g=1.164*(s1.val[0]-16)-0.813*(s3.val[0]-128)-0.391*(s2.val[0] - 128);
		r=1.164*(s1.val[0]-16) + 1.596*(s3.val[0]-128);
		s.val[0]=b;
		s.val[1]=g;
		s.val[2]=r;
		cvSet2D(dst,i,j,s);
	}
cvSaveImage("./ResultImg/result.bmp",dst);
}


//Apply Log to origional image
void selflog(IplImage *src, IplImage *dst)
{
IplImage *src1=0;
src1 = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
cvConvertScale(src, src1, 1, 0);
cvAddS( src1, cvScalar(1.0), src1 );
cvLog( src1,dst ); 
}


//Apply Exp at the last step 
void selfexp(IplImage *src, IplImage *dst)
{
cvExp( src,dst );
cvAddS( dst, cvScalar(-1.0), dst );
}


//Forward Fourier Trans
void fft2(IplImage *src, IplImage *dst)
{
int height = src->height;
int width = src->width;
IplImage *src1=0,*image_Re = 0, *image_Im = 0, *Fourier = 0;
image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); 
image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); 
Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);
src1 = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
cvConvertScale(src, src1, 1, 0);
for(int j=0;j<width;j++){
	  for(int i=0;i<height;i++){
		  int mod=(i+j)%2;
		  if (mod==0)
		  CV_IMAGE_ELEM(src1,double,j,i)=CV_IMAGE_ELEM(src1,double,j,i);
		  else
		  CV_IMAGE_ELEM(src1,double,j,i)=(-1)*CV_IMAGE_ELEM(src1,double,j,i);
	  }
	 }

cvConvertScale(src1, image_Re, 1, 0);
cvZero(image_Im);
cvMerge(image_Re, image_Im, 0, 0, Fourier);

// Application of the forward Fourier transform
cvDFT(Fourier, dst, CV_DXT_FORWARD);
cvReleaseImage(&image_Re);
cvReleaseImage(&image_Im);
cvReleaseImage(&Fourier);
}


void fft2shift(IplImage *src, IplImage *dst)
{
IplImage *image_Re = 0, *image_Im = 0;
int nRow, nCol, i, j, cy, cx;
double scale, shift;
image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
cvSplit( src, image_Re, image_Im, 0, 0 );

cvPow( image_Re, image_Re, 2.0);
cvPow( image_Im, image_Im, 2.0);
cvAdd( image_Re, image_Im, image_Re);
cvPow( image_Re, image_Re, 0.5 );

cvAddS( image_Re, cvScalar(1.0), image_Re ); 
cvLog( image_Re, image_Re ); 

// Normalize image (0 - 255)
double minVal = 0, maxVal = 0;
cvMinMaxLoc( image_Re, &minVal, &maxVal );
scale = 255/(maxVal - minVal);
shift = -minVal * scale;
cvConvertScale(image_Re, dst, scale, shift);
cvReleaseImage(&image_Re);
cvReleaseImage(&image_Im);
}


int main()
{
IplImage *srcCol; 
IplImage *src; 
IplImage *srcmid;
IplImage *Fourier;
IplImage *dst ;
IplImage *ImageRe;
IplImage *ImageIm;
IplImage *ImageRe1;
IplImage *ImageIm1;
IplImage *Image;
IplImage *ImageDst;

srcCol = cvLoadImage("./TestImg/Dom.png",1);
IplImage* imgYUV=cvCreateImage(cvGetSize(srcCol), IPL_DEPTH_64F, 3);
RGB2YUV(srcCol, imgYUV);

double m,M;
double scale;
double shift;
//src = cvLoadImage("./TestImg/Dom.png",0);
src = cvLoadImage("./TestImg/out_y.bmp",0);
cvNamedWindow("Source",1);
cvShowImage("Source",srcCol);             

int height = src->height;
int width = src->width;

//Log
srcmid = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
selflog(src,srcmid);

Fourier = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
dst = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
ImageRe = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
ImageIm = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
ImageRe1 = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
ImageIm1 = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
Image = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
ImageDst = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);

//Fourier Trans
fft2(srcmid,Fourier);
fft2shift(Fourier, Image); 
cvSplit(Fourier,ImageRe1,ImageIm1,0,0);


//Impliment Homomorphic filters
double rl=0.2;
double rh=3;
double c=0.2;

for(int j=0;j<width;j++){
	  for(int i=0;i<height;i++){
		  double mid1=1-exp(c*((i-height/2)*(i-height/2)+(j-width/2)*(j-width/2))/((0-height/2)*(0-height/2)+(0-width/2)*(0-width/2)));
		  double mid=(rh-rl)*mid1+rl;
		  CV_IMAGE_ELEM(ImageRe1,double,j,i)=mid*(CV_IMAGE_ELEM(ImageRe1,double,j,i));
		  CV_IMAGE_ELEM(ImageIm1,double,j,i)=mid*(CV_IMAGE_ELEM(ImageIm1,double,j,i));
		  }
}

cvMerge(ImageRe1, ImageIm1, 0, 0, Fourier);

//Inverse Fourier Trans
cvDFT(Fourier,dst,CV_DXT_INV_SCALE);
cvSplit(dst,ImageRe,ImageIm,0,0);

for(int j=0;j<width;j++){
	  for(int i=0;i<height;i++){
		  int mod=(i+j)%2;
		  if (mod==0)
		  CV_IMAGE_ELEM(ImageRe,double,j,i)=CV_IMAGE_ELEM(ImageRe,double,j,i);
		  else
		  CV_IMAGE_ELEM(ImageRe,double,j,i)=(-1)*CV_IMAGE_ELEM(ImageRe,double,j,i);
	  }
	 }
//cvPow(ImageRe,ImageRe,2);               
//cvPow(ImageIm,ImageIm,2);
//cvAdd(ImageRe,ImageIm,ImageRe,NULL);
//cvPow(ImageRe,ImageRe,0.5);

//Exp
selfexp( ImageRe, ImageRe );

//Normalization
cvMinMaxLoc(ImageRe,&m,&M,NULL,NULL);
scale = 235/(M - m);
shift = -m * scale;
cvConvertScale(ImageRe,ImageDst,scale,shift);

//cvNamedWindow("Fourier",0);
//cvShowImage("Fourier",Image);
//cvNamedWindow("Inversed",0);
//cvShowImage("Inversed",ImageDst);

//cvSaveImage("./ResultImg/Fourier.bmp",Image);
//cvSaveImage("./ResultImg/Dst.bmp",ImageDst);
cvSaveImage("./TestImg/after_y.bmp",ImageDst);

IplImage *src1,*src2,*src3; 
src1 = cvLoadImage("./TestImg/after_y.bmp",1);
src2 = cvLoadImage("./TestImg/out_u.bmp",1);
src3 = cvLoadImage("./TestImg/out_v.bmp",1);
IplImage* imgRGB=cvCreateImage(cvGetSize(src1), IPL_DEPTH_8U, 3);
YUV2RGB(src1, src2, src3, imgRGB);

cvNamedWindow("Inversed",3);
cvShowImage("Inversed",imgRGB);

//Release
cvWaitKey(1000000);
cvReleaseImage(&src);
cvReleaseImage(&Image);
cvReleaseImage(&ImageIm);
cvReleaseImage(&ImageRe);
cvReleaseImage(&Fourier);
cvReleaseImage(&dst);
cvReleaseImage(&ImageDst);
return 0;
}

